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Abstract We consider the theory of fluctuations of a colloidal solid observed in a confocal slice. For a 
cubic crystal we study the evolution of the projected elastic properties as a function of the anisotropy 
of the crystal using numerical methods based on the fast Fourier transform. In certain situations of high 
symmetry we find exact analytic results for the projected fluctuations. 
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1 Introduction 

In a recent paper [T] we simulated a large ensemble of hard 
spheres packed in a face-centered cubic crystal in order to 
study the fluctuations within a slice of the sample. This 
work was motivated by the interpretation of the elastic 
properties found in experiments of colloidal crystals OE1 
4,5,6,7 . In particular we showed how to extract informa- 
tion on the three-dimensional elastic constants from the 
anomalous dispersion relationship that is observed within 
a single slice. While we managed to perform an exact an- 
alytic calculation for isotropic elastic media many exper- 
iments are performed on crystals, rather than amorphous 
isotropic solids. The question then arises as to how to best 
treat the anisotropy. We proposed using an averaging pro- 
cedure such as that of Fedorov [Hl[9] and for the sample 
that we simulated the results were remarkably good: Simu- 
lation and theory agreed to within a few percent. However, 
many months of simulation are required to study a single 
sample and it is difficult to perform a systematic study of 
the effect of anisotropy on the projected properties. 

In this paper we study the projected fluctuations within 
linear elasticity, rather than by molecular dynamics sim- 
ulation. By formulating cubic elastic properties on a dis- 
crete grid we are able to evaluate the anisotropic Green 
function within a projection plane using fast Fourier trans- 
forms. This procedure allows us to study how well our an- 
alytic approximation reproduces the exact projected prop- 
erties as a function of the plane orientation as well as of 
the three elastic constants characterizing the full cubic 
system. 



2 Cubic Elasticity 

Let us now consider a three-dimensional cubic crystal of 
atoms, indexed by i, deformed by the displacement vec- 
tor Ui. We write the elastic energy as a quadratic form in 



the spatial derivatives utj of this vector. This quadratic 
form, which respects the symmetry of the crystal, can be 
expressed with the help of the Christoffel matrix [TU] 

Di k (k) = A6 i j5 kl + n(5 ik 5j l + 6 u 5 jk )+vSi ikl kjk : , (1) 

with a reciprocal vector k, Lame constants A, |x, and an- 
isotropy v. The tensor S = Y.-p=i e p e p e p e p , with e p unit 
vectors parallel to the cubic axes of the crystal. Summa- 
tion over indices occurring twice is assumed. In a colloidal 
crystal, being under external pressure P, the constants are 
modified by non-linear interactions [TO]. The above con- 
stants are then related to the elastic constants used in 
Voigt notation, 

C u -P=A + 2|i + -v 

C44 - P = M- 
Ci 2 +P = A. 

The Green function of the static elastic problem is then 
the inverse of the Christoffel matrix, 



D ij (k)G jlc (k) = & ij . 



(2) 



One expresses the free energy in terms of the displacement 
field 

F[u] = ^u i (k)D lj (k)u j (-k). (3) 



For each wavevector k, D is a three-by-three matrix with 
eigenvalues dt(k) where the subscript i indicates a polar- 
ization state. Following a convention usual in the liter- 
ature [miSlIT]. we define the auxiliary variable tu?(k) = 
di'k:. 

In our previous work we showed that it was useful to 
perform an angular average of the elastic constants in or- 
der to generate the "best" isotropic approximation to the 
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elastic properties [HE]. This average gives the effective 
Lame constants 

~v ~v 

A = A + - and p. = [L + - . (4) 

5 5 

On observing a slice of the three-dimensional solid un- 
der a confocal microscope one can then show that the 
effective dispersion is described by an anomalous rela- 
tion where iv ~ q 1 / 2 , for small q pQ. Here q is the two- 
dimensional wavevector, and is to be distinguished from 
the full three-dimensional wavevector k. A detailed calcu- 
lation gives 

iu\ — 2Q.q (transverse), 

2 4a(A+2u) . (5) 

cum = — q (longitudinal). 

A + 3p. 

Thus observation of the two branches of the projected 
dispersion curve gives information on the averaged Lame 
moduli. It is interesting to note that even in the incom- 
pressible limit of rubber elasticity where u/A — > the ef- 
fective projected frequencies remain finite. Apparent com- 
pression in the two-dimensional plane remains easy. In this 
limit we find that cu^/cu 2 = 1/2. 

3 Discretization 

Rather than performing a molecular dynamics simulation 
of hard spheres, we here use only methods from linear al- 
gebra to numerically project and characterize fluctuations 
from an elastic solid. The disadvantage is that all non- 
linearities in the true physical system are neglected. The 
advantage is that very efficient methods are available that 
generate useful data in just a few minutes of calculation. 
One is also able to systematically vary the parameters of 
the model in order to study their influence on the pro- 
jected fluctuations. 

We use a finite difference discretization of Eq. ([I]), 
taking the mesh size as unity and embedding the elastic 
medium in a cubic box of dimensions L. For any wavevec- 
tor k we define the components of a vector v as 

Vj=(l-e ik a (6) 

The finite different discretization corresponding to Eq. ([I]) 
is then 
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Figure 1. Longitudinal and transverse fluctuations of 
an isotropic system projected to the two-dimensional 
plane (1,0,0). The result is compared to the long-wavelength 
limit (dotted line) predicted in Eq. L = 512, A = 70, u = 95, 
v = 0. 

are found by rotating the crystal with respect to the cu- 
bic cell, modifying the contribution in v. This is easily 
performed by updating the tensor S using rotation ma- 
trices parameterized with Euler angles. In this case the 
vectors e p are the three rows of the Euler rotation matrix 
and we construct the anisotropic elastic contribution from 
|eP.v| 2 e?e?. 

We evaluate D(k) for wavevectors ki. = 27Trii/L. We 
numerically invert the three-by-three matrix D(k) for each 
value of k and then use the fast Fourier transform to eval- 
uate the components of the real space Green function on a 
cubic grid. We then evaluate the Green functions on a slice 
and use two-dimensional fast Fourier transforms to obtain 
the effective dispersion relation for a projected wavevec- 
tor q. Such a dispersion relation is shown in Fig. [T] for the 
isotropic case v = 0, where we compare with the contin- 
uum expressions Eq. ([5]). As expected, we see agreement 
between theory and simulation. We use values of u. and A 
comparable to those found in our previous molecular dy- 
namics simulations. The anisotropy v will be varied over 
a wide range of values in the figures below. We use units 
such that the energy scale is k^T. 

We now show how to calculate the full dispersion curve 
corresponding to the present discretization, using a alter- 
native calculational route than that proposed in our orig- 
inal paper. 



D = uz(k) 1+ (A+n.)(v(8> v*) +v diag(vivt), (7) 

withz(k) = tr(v<g>v*) = ^ t (2— 2 cos k^) and where "diag" 
constructs a three-by-three matrix with the given compo- 
nents on the diagonal. I is the three-dimensional unit ma- 
trix, and ® denotes the exterior product of two vectors. 
The matrix D is complex Hermitian. One sees that the 
small-wavevector expansions of Eq. Q is indeed given by 
the Christoffel matrix. The exact analytic form makes the 
dispersion relation periodic within the first Brillouin zone 
of the discretizing mesh. Note this form of the discretiza- 
tion suffices to study cuts in the plane (1, 0, 0). Other cuts 



3.1 Projection for v = 

In our previous work we gave a treatment of the projec- 
tion which required calculation of the explicit form of the 
real-space Green function. We show here that a different 
method allows one to calculate the projection for v = 0, 
and we apply this method to the discretized Eq. ^ , gen- 
erating the analytic expression for the full curves of Fig. [I] 
We firstly calculate the inverse of Eq. 

G(k) = ^(k)v -A + 2u zfk) J' (8) 
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We now formally calculate the real-space Green function 
by Fourier transformation, and back transform the pro- 
jected correlations. Elementary calculations show that 



G 2 (q) = ^G(q x ,q y ,k z 



1 

2tt 



dk z G(q, k z 



(9) 



where G2 describes the correlations projected on the con- 
focal plane z = 0. 

The integral is calculated by noting that 



1 

2n 



1 



1 

27T 



A — 2 cos k z 
1 



-dk z 



1 



_ n (A — 2 cosk z 



rdk z = 



(A 2 -4)V2- 
A 

(A 2 - 4) 3 / 2 ' 



Within the plane we consider without loss of generality 
the choice q — ( q , 0) , then A = 4 — 2 cos q . The dispersion 
relations are then 



1 



1 



1 



CO 



_l H [(4 -2 cos q) 2 -4] 



1/2 



(10) 



1 1 A+p. / (4-2cosq)(2-2cosq) 



con 



liA + 2u ^ [(4 _ 2 cos q) 2 -4] 



3/2 



When one expands Eq. ( 10 1 about q = one finds exactly 
Eq. It should be noted that all higher order correc- 
tions in the expansion of Eq. ( 10 ) are non-universal, and 
thus change on using different discretizations of the con- 
tinuum equations. However these equations are useful in 
calibrating the convergence of the numerical code that we 
developed for this paper. 



box, for which L z = L the image interaction via this com- 
ponent of the density is strongly suppressed. Higher order 
waves have even weaker interactions. 

We now return to the problem of an isotropic elastic 
solid. If we observe a sinusoidal perturbation on the plane 
z = we can find the minimum energy fluctuation con- 
sistent with this observation by examining the solution of 
the equation 



(1 — 2ct)V"u + graddivu = 



(12) 



away from the plane z = 0, where cr is the Poisson ratio 
of the material. We look for solutions of the form 



-ikx 1 



u(z),0,w(z)) 



(13) 



and find 



(1 - 2a) (-k 2 + y 2 )u - k 2 u + ikyw = 0, 
(1 - 2cr)(-k 2 +y 2 )w + y 2 w + ikyu = 0, 

where y is the Laplace transform of the fields in the z 
direction. Solutions are found when the determinant of 
coefficients is zero which gives 

k 2 =y 2 , 

implying a result very similar to that in electrostatics: 
There is a decay of interactions on a scale which is de- 
termined by the wavelength of the observed mode, inde- 
pendent of the Poisson ratio. There is also a simple phase 
relationship between the horizontal and vertical compo- 
nents of the fluctuations: 



4 Three-dimensional structure of fluctuations 

In this section we discuss the three-dimensional nature of 
the fluctuation field induced in the presence of a fluctu- 
ation in a surface. In particular we wish to argue that 
the use of periodic boundary conditions does not lead to 
strong artifacts in the amplitudes of the measured waves, 
at the same time we are able to estimate the decay of 
boundary perturbations to a sample on fluctuations within 
a sample. Very similar arguments are regularly used in 
electrostatics [T2"lll3|. Let use firstly summarize the sit- 
uation for solutions of the Laplace equation where the 
algebra is particularly simple. 

Consider a periodic box of lateral dimensions L and 
height L z . Let there be charges on the plane z = mod- 
ulated with wavevector k x . Away from the plane, in the 
absence of sources, the potential is given by a solution to 
the equation V 2 cf) = 0. If the solution is periodic in the 
x — y plane then the solution must decay exponentially in 
the z direction so that 

4> = cf) k exp (ilex — k|z|). (11) 

The slowest decaying wave in the z direction has k = 27t/L, 
thus successive images in the box have an asymptotic in- 
teraction [12 that varies as ~ e~ 27tLz / L . Already in a cubic 



u ± iw = 

From this calculation we understand that we should 
only observe samples far from any external walls. Fluctu- 
ations of wavelength i observed within a plane decay over 
a distance l/2n and will be perturbed by external surfaces 
which are too close. This result allows us to find a simple 
physical interpretation of the anomalous dispersion rela- 
tion: Imposition of a fluctuation of wavelength I leads to 
an energy density (u/£) 2 . This excites a total volume of 
the sample which varies as L 2 £, giving a scaling in the total 
elastic energy of mode q as F q ~ 1/1 ~ | q | . 

We now study the variation of the elastic properties 
as a function of the cubic anisotropy v, and compare the 
results with the averaging approximation. 



5 Numerical Results 

We implemented the Fourier and matrix manipulation us- 
ing an Octave/Matlab script; the language was chosen for 
its facilities for the manipulation of fast Fourier transforms 
and matrix algebra. We projected the fluctuations onto 
different planes of the crystal and studied the effective 
dispersion relations, Fig. [2j For the cut of the crystal in 
the plane (1, 0, 0) we plot the effective dispersion relations 
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Figure 2. Evolution of amplitudes as a function of y in a 
cut perpendicular to (1,0,0). Previous theory in dotted line, 
Eq. (J5|. The theoretical curves are intermediate between the 
two numerically generated curves. The curve in the direction 
(1,0) displays the exceptional features to be independent of v. 
All lines intersect at v = since the system becomes isotropic. 
L = 512. A/^ = 0.73. 



useful experimentally, because the samples naturally grow 
from fiat surfaces aligned in this orientation. 
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Figure 4. Evolution of amplitudes as a function of "V in a cut 
perpendicular to (1,1,1). The dispersion curve is isotropic in 
the plane. L = 512. A/u. = 0.73. 



in two directions in the plane, (1,0) and (1,1). The av- 
eraging approximation, while not wholly inaccurate, does 
not account for a number of qualitative features of the 
measured dispersion curves. In particular, the measured 
transverse stiffness in the direction (1,1) is independent 
of v while the longitudinal mode varies strongly with the 
anisotropy. This gives rise to level crossing for large nega- 
tive values of v, which occurs for values of the parameters 
very close to those found for the hard sphere system. 

The situation in the plane (1,1,0), Fig.[3j is somewhat 
closer to the approximate analytic form ofthe dispersion 
curves. Here again the dispersion is anisotropic within the 
plane. Clearly this is impossible to describe with spheri- 
cally averaged elastic constants. 



The vertical gray lines in Figs.[2j|4|at v/u. = — 1.06 rep- 
resent the anisotropy found in the full molecular-dynamics 
simulation. Also the ratio has been chosen to fit to the 
simulation. Comparing the published data pQ with Fig. [4j 
we notice that the full simulation curves are displaced with 
respect to the data produced by Fourier analysis. Both are 
slightly displaced with respect to the isotropic average, 
but with different signs. We interpreted this shift as being 
due to non-linearities in the full system and estimated it to 
about 3%. A new estimate based on the numerical method 
described in this paper is larger, about 9% and 15%. 



5.1 Analytic projection on (1,0,0) 
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Figure 3. Evolution of amplitudes as a function of v in a cut 
perpendicular to (1, 1, 0). L = 512. A/u. = 0.73. 



Using the method of Sec. |3j it is easy to write a for- 
mal solution for the projected properties for arbitrary y. 
In the limit of small q one finds as integrand the ra- 
tio of two polynomials of orders four and six in the vec- 
tor (q x ,q y ,k z ), see also [H]. In principle by factorizing 
the denominator (which is of order six) one has the full 
analytic solution to the problem. In practice this requires 
the roots of high order equations and leads to clumsy and 
opaque expressions. There is, however, one case which is 
particularly easy to understand, namely the independence 
of the transverse mode in Fig. [3] on y. This is a sim- 
ple consequence of the Sherman-Morrison identity: The 
Christoffel matrix on the plane (1, 0, 0) can be written as 
a diagonal matrix plus a rank one update: 



with 



D = B + (A + |x)(k<8k), 
B = diag(uk 2 +vkf). 



For the plane (1, 1, 1), Fig. [3]we see that the averaging 
procedure works particularly well and that the dispersion 
relations are isotropic. This is the curve which is the most 



The inverse of this matrix is then 



(A 



+ uO(B- 1 k)®(B- 1 k) 
1 + (A + u.)kB-!k 



(14) 
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If q = (1, 0), then it is clear that G xy (1, 0, k z ) = and the 
only contribution to the projected transverse spectrum 
comes from G y y (1, 0, k z ) which is trivially independent 
of v. 

The longitudinal contribution can also be found from 
the Sherman-Morrison formula, and we find the integral: 



dk 



1 



_ n 2tt1 + v/li + k 2 



1 + k 2 + k 2 v/u- 



(l+v/|x+k 2 )(l+k 2 +k 2 v/nQ 
1+A/vi 



+ l + 2(l + v/u)k 2 + k 4 



This can be integrated by finding the roots of the denom- 
inator, and expanding using partial fractions, the result is 
however too complicated to be illuminating. 



5.2 Simulation data 

On seeing the surprising nature of the dispersion in Fig. [2] 
where the longitudinal and transverse modes become al- 
most degenerate at large negative values of the anisotropy, 

around v (X, we performed an analysis based on our 

molecular dynamics simulation, Fig. [5] There, the aniso- 
tropy was also strong, v = — 1.06u. We find that indeed 
our numerical data displays very similar features, confirm- 
ing the useful nature of the Fourier based code for explor- 
ing the properties of the elastic behavior of hard sphere 
systems. The longitudinal and transverse dispersion rela- 
tions are almost degenerate over all wavevectors in the 
direction (1,0). We do not have any simple explanation 
for this result. 



2() 



10 



3" 5 



2 - 




longitudinal 
transverse 
14.79^fk| 
12.24 Y/jkf 



0.05 0.1 



0.2 



0.5 
|k| 



10 



Figure 5. Effective dispersion relations for the plane (1,0,0) 
in the direction (1,0). Dataset of pQ. The longitudinal and 
transverse branches are almost degenerate for large negative 
cubic anisotropy. 



6 Density of states 

In the interpretation of the density of states of disordered 
material one often compares with the Debye theory of 




Figure 6. Density of states measured on the (1,1,1) cut, 
data from the molecular dynamics simulation Q] . (a) Note the 
prominent van Hove singularities corresponding to the longi- 
tudinal and transverse projected modes. Panel (b) plots the 
integrated density of states on a log-log scale, compared with 
the theoretical form cu 4 . 



elasticity. This allows one to establish an expectation for 
the low-frequency mode structure, and then to determine 
whether the disorder has given rise to a deficit or excess 
in the density of states. As noted in [15] a dispersion re- 
lation such as Eq. ^ leads to modifications in the Debye 
spectrum, which is normally given as p(cu) - cu d_1 for a 
d-dimensional elastic solid. If we use instead the dispersion 
relations in Eq. ([5| and define cu 2 = a^q and w\ — a±q 
we find that the density of states, expressed in terms of 
the variable cu is given by 



p(cu) = 



1 / 1 



cu' 



(15) 



which corresponds to neither the projected, nor the em- 
bedding dimension of the slice. Figure [6] shows the density 
of states of a cut of our molecular dynamics crystal. The 
top panel displays the binned density of states on a linear 
scale. In order to avoid artifacts of binning we also plot 
the integrated density of states in the bottom panel and 
find that there is a good fit to a law in cu 4 for cu < 8. Some 
of the experimental literature plots two-dimensional slices 
normalized by the Debye density of states for a three- 
dimensional crystal [5]. We have also plotted this form 
in the top panel of Fig. [7j We see, rather surprisingly, a 
region between the local frequency behavior in cu 3 and 
the first peak where the normalized data seems relatively 
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Figure 7. Density of states divided by cu 2 , projected on two 
different planes. The resulting curves depend on the projection 
plane: In the case (1,1,1) the expected regime in p - cu 3 is seen 
for cu < 10; a plateau is observed for 10 < cu < 20. The form 
p ~ cu 3 fits over a larger range of cu for the case (1, 0, 0). 



constant. This seems, however, to be a simple cross-over 
occurring over a modest range of frequencies. If we project 
on a different plane, (1, 0, 0) in the bottom panel of Fig.JTJ 
no plateau is visible. 
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7 Conclusions 

We have introduced a Fourier based method which can ef- 
ficiently generate the Green function in anisotropic elastic 
media. It allows one to rapidly perform detailed studies of 
the evolution of correlations and fluctuations in projected 
geometries, such as those used in several recent experi- 
ments. We have studied the effective dispersion relations 
of a colloidal crystal projected on the plane (1,0,0) and 
shown that the near degeneracy of the dispersion curves in 
the direction (1,0) is due to the large negative cubic aniso- 
tropy. Surprisingly long wavelengths are needed to see the 
continuum limit in the density of states. Studies on cer- 
tain cuts lead to an intermediate regime for the density of 
states. 

Code written in Octave/Matlab is available as sup- 
plementary material to the paper to perform the Fourier 
based evaluation of Green functions, together with the 
projection to two dimensions. 
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